Additional Parameters
ISI
# Test ISI
dat <- data.frame()
for(i in c("Ebbinghaus", "MullerLyer", "VerticalHorizontal")) {
dat <- rbind(
mgcv::gamm(RT ~ s(ISI),
random = list(Participant = ~1),
data=filter(perceptual, Illusion_Type == i, Error==0)) |>
modelbased::estimate_relation(length=30) |>
mutate(Illusion_Type = i, Outcome = "RT", type="GAM"),
glmmTMB::glmmTMB(RT ~ poly(ISI, 2) + (1|Participant),
data =filter(perceptual, Illusion_Type == i, Error==0)) |>
modelbased::estimate_relation(length=30) |>
select(-Participant) |>
mutate(Illusion_Type = i, Outcome = "RT", type="poly"),
mgcv::gamm(Error ~ s(ISI),
random = list(Participant = ~1),
data=filter(perceptual, Illusion_Type == i),
family = "binomial") |>
modelbased::estimate_relation(length=30) |>
mutate(Illusion_Type = i, Outcome = "Error", type="GAM"),
glmmTMB::glmmTMB(Error ~ poly(ISI, 2) + (1|Participant),
data =filter(perceptual, Illusion_Type == i),
family = "binomial") |>
modelbased::estimate_relation(length=30) |>
select(-Participant) |>
mutate(Illusion_Type = i, Outcome = "Error", type="poly")
) |>
rbind(dat)
}
dat |>
ggplot(aes(y = Predicted, x=ISI)) +
geom_ribbon(aes(ymin=CI_low, ymax=CI_high, fill = type), alpha=0.3) +
geom_line(aes(color = type)) +
facet_wrap(Outcome ~ Illusion_Type, scales = "free")
Ebbinghaus
Model Selection
test_models <- function(data) {
# TODO: add random effect
models_err <- list()
models_rt <- list()
for(f in c("Illusion_Difference",
"logmod(Illusion_Difference)",
"sqrtmod(Illusion_Difference)",
"cbrtmod(Illusion_Difference)")) {
err <- glmmTMB::glmmTMB(as.formula(paste0("Error ~ ", f, "+ (1|Participant)")),
data = data, family = "binomial")
models_err[[f]] <- err
rt <- glmmTMB::glmmTMB(as.formula(paste0("RT ~ ", f, " + poly(ISI, 2) + (1|Participant)")),
data = filter(data, Error == 0))
models_rt[[f]] <- rt
}
mutate(performance::test_performance(models_err), Outcome = "Error") |>
rbind(mutate(performance::test_performance(models_rt), Outcome = "RT")) |>
select(-Model, -log_BF) |>
datawizard::convert_na_to(select="BF", replacement = 1) |>
arrange(Outcome, desc(BF)) |>
export_table(footer = "Each model is compared to 'Illusion_Difference'")
}
test_models(filter(perceptual, Illusion_Type == "Ebbinghaus"))
Error Rate
data <- filter(perceptual, Illusion_Type == "Ebbinghaus")
Descriptive
plot_desc_errors <- function(data) {
data |>
ggplot(aes(x = Illusion_Difference)) +
geom_histogram(data=filter(data, Error == 1),
aes(y=..count../sum(..count..), fill = Illusion_Side),
binwidth = diff(range(data$Illusion_Difference)) / 20, color = "white") +
geom_smooth(aes(y = Error, color = Illusion_Side),
method = 'gam',
formula = y ~ s(x, bs = "cs"),
method.args = list(family = "binomial")) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
scale_color_manual(values = c("-1" = "#FF5722", "1" = "#43A047")) +
scale_fill_manual(values = c("-1" = "#FF5722", "1" = "#43A047")) +
coord_cartesian(ylim = c(0, 1), xlim = range(data$Illusion_Difference)) +
labs(x = "Task Difficulty", y = "Probability of Error") +
theme_modern()
}
plot_desc_errors(data)

Model Specification
formula <- brms::bf(
Error ~ sqrtmod(Illusion_Difference) +
(1 + sqrtmod(Illusion_Difference) | Participant),
family = "bernoulli",
decomp = "QR"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
perceptual_ebbinghaus_err <- brms::brm(formula,
data = data,
refresh = 0,
normalize = FALSE
)
save(perceptual_ebbinghaus_err, file="models/perceptual_ebbinghaus_err.Rdata")
Model Inspection
load("models/perceptual_ebbinghaus_err.Rdata")
plot_model_errors <- function(data, model) {
pred <- estimate_relation(model, at="Illusion_Difference", length = 100)
data |>
ggplot(aes(x = Illusion_Difference)) +
geom_histogram(data=filter(data, Error == 1),
aes(y=..count../sum(..count..)),
binwidth = diff(range(data$Illusion_Difference)) / 20) +
geom_ribbon(data = pred,
aes(ymin = CI_low, ymax = CI_high),
alpha = 1/3, fill = "red") +
geom_line(data = pred,
aes(y = Predicted),
color = "red") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
coord_cartesian(ylim = c(0, 1), xlim = range(data$Illusion_Difference)) +
labs(x = "Task Difficulty", y = "Probability of Error") +
theme_modern()
}
plot_model_errors(data, perceptual_ebbinghaus_err)

Response Time
data <- filter(perceptual, Illusion_Type == "Ebbinghaus", Error == 0)
Descriptive
plot_desc_rt <- function(data) {
data |>
ggplot(aes(x = Illusion_Difference, y = RT)) +
# ggpointdensity::geom_pointdensity(size = 3, alpha=0.5) +
# scale_color_gradientn(colors = c("grey", "black"), guide = "none") +
# ggnewscale::new_scale_color() +
stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
scale_fill_gradientn(colors = c("white", "black"), guide = "none") +
ggnewscale::new_scale_fill() +
geom_smooth(aes(color = Illusion_Side, fill = Illusion_Side),
method = 'gam',
formula = y ~ s(x, bs = "cs")) +
scale_color_manual(values = c("-1" = "#FF5722", "1" = "#43A047")) +
scale_fill_manual(values = c("-1" = "#FF5722", "1" = "#43A047")) +
scale_x_discrete(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = "Task Difficulty", y = "Response Time (s)") +
coord_cartesian(ylim = c(0, 2.5)) +
theme_modern() +
ggside::geom_ysidedensity(aes(fill = Illusion_Side), color = NA, alpha = 0.3) +
ggside::theme_ggside_void() +
ggside::scale_ysidex_continuous(expand = c(0, 0)) +
ggside::ggside()
}
plot_desc_rt(data)

Model Specification
# TODO: Add random to parameters
formula <- brms::bf(
RT ~ sqrtmod(Illusion_Difference) + poly(ISI, 2) +
(1 + sqrtmod(Illusion_Difference)| Participant),
sigma ~ sqrtmod(Illusion_Difference) + poly(ISI, 2) +
(1 + sqrtmod(Illusion_Difference)| Participant),
beta ~ sqrtmod(Illusion_Difference) + poly(ISI, 2) +
(1 + sqrtmod(Illusion_Difference)| Participant),
family = "exgaussian",
decomp = "QR"
)
# brms::get_prior(formula, data = data)
# brms::validate_prior(formula)
perceptual_ebbinghaus_rt <- brms::brm(formula,
data = data,
refresh = 0,
init = 0,
normalize = FALSE
)
save(perceptual_ebbinghaus_rt, file="models/perceptual_ebbinghaus_rt.Rdata")
Model Inspection
load("models/perceptual_ebbinghaus_rt.Rdata")
plot_model_rt <- function(data, model) {
pred <- estimate_relation(model, at="Illusion_Difference", length = 100)
data |>
ggplot(aes(x = Illusion_Difference)) +
stat_density_2d(aes(fill = ..density.., y = RT), geom = "raster", contour = FALSE) +
scale_fill_gradientn(colors = c("white", "black"), guide = "none") +
ggnewscale::new_scale_fill() +
geom_ribbon(data = pred,
aes(ymin = CI_low, ymax = CI_high),
alpha = 1/3, fill = "red") +
geom_line(data = pred,
aes(y = Predicted),
color = "red") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = "Task Difficulty", y = "Response Time (s)") +
coord_cartesian(ylim = c(0, 2.5)) +
theme_modern()
}
plot_model_rt(data, perceptual_ebbinghaus_rt)

Individual Scores
get_scores <- function(model, illusion="Ebbinghaus") {
family <- insight::find_response(model)
scores <- modelbased::estimate_grouplevel(model) |>
data_filter(str_detect(Level, "Participant")) |>
mutate(Group = str_remove(Group, ": Participant"),
Level = str_remove(Level, "Participant."),
Level = str_remove(Level, "_beta."),
Level = str_remove(Level, "_sigma."),
Group = str_remove_all(Group, "cbrtmod"),
Group = str_remove_all(Group, "sqrtmod"),
Group = str_remove_all(Group, "sqrt"),
Group = str_remove_all(Group, "logmod"),
Group = str_remove_all(Group, "abs"),
Group = str_replace(Group, "Illusion_Difference", "Diff"),
Group = str_replace(Group, "Intercept__sigma", "InterceptSigma"),
Group = str_replace(Group, "Intercept__beta", "InterceptBeta"),
Group = str_replace(Group, "Diff__sigma", "DiffSigma"),
Group = str_replace(Group, "Diff__beta", "DiffBeta"))
p <- scores |>
ggplot(aes(x = Median, y = Level)) +
geom_pointrange(aes(xmin = CI_low, xmax = CI_high, color = Group), size=0.5) +
scale_color_flat_d(guide = "none") +
scale_fill_flat_d(guide = "none") +
labs(y = "Participants") +
theme_modern() +
theme(strip.placement = "oustide",
axis.title.x = element_blank(),
axis.text.y = element_blank()) +
ggside::geom_xsidedensity(aes(fill=Group, y = after_stat(scaled)), color = NA, alpha = 0.3) +
ggside::theme_ggside_void() +
ggside::scale_xsidey_continuous(expand = c(0, 0)) +
ggside::ggside() +
facet_grid(~Group, switch = "both", scales = "free") +
ggtitle(paste(illusion, "-", family))
scores <- scores |>
select(Group, Participant = Level, Median) |>
pivot_wider(names_from = "Group", values_from = "Median") |>
data_rename("Diff",
paste0("Perception_", illusion, "_Difficulty_", family)) |>
data_rename("DiffSigma",
paste0("Perception_", illusion, "_Difficulty_SigmaRT")) |>
data_rename("DiffBeta",
paste0("Perception_", illusion, "_Difficulty_BetaRT")) |>
data_rename("Intercept",
paste0("Perception_", illusion, "_Intercept_", family)) |>
data_rename("InterceptSigma",
paste0("Perception_", illusion, "_Intercept_SigmaRT")) |>
data_rename("InterceptBeta",
paste0("Perception_", illusion, "_Intercept_BetaRT"))
list(scores = scores, p = p)
}
ebbinghaus_err <- get_scores(perceptual_ebbinghaus_err, illusion="Ebbinghaus")
ebbinghaus_rt <- get_scores(perceptual_ebbinghaus_rt, illusion="Ebbinghaus")
mullerlyer_err <- get_scores(perceptual_mullerlyer_err, illusion="MullerLyer")
mullerlyer_rt <- get_scores(perceptual_mullerlyer_rt, illusion="MullerLyer")
verticalhorizontal_err <- get_scores(perceptual_verticalhorizontal_err, illusion="VerticalHorizontal")
verticalhorizontal_rt <- get_scores(perceptual_verticalhorizontal_rt, illusion="VerticalHorizontal")
p <- (ebbinghaus_err$p + ebbinghaus_rt$p) /
(mullerlyer_err$p + mullerlyer_rt$p) /
(verticalhorizontal_err$p + verticalhorizontal_rt$p) +
plot_layout(guides = "collect") +
plot_annotation(title = "Inter- and Intra- Variability of Perceptual Scores", theme = theme(plot.title = element_text(face = "bold", hjust = 0.5)))
p

scores <- ebbinghaus_err$scores |>
merge(ebbinghaus_rt$scores, by="Participant") |>
merge(mullerlyer_err$scores, by="Participant") |>
merge(mullerlyer_rt$scores, by="Participant") |>
merge(verticalhorizontal_err$scores, by="Participant") |>
merge(verticalhorizontal_rt$scores, by="Participant")
write.csv(scores, "../data/scores_perceptual.csv", row.names = FALSE)